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Dynamical properties of the impurity spin-| in 2D and quasi-2D Heisenberg antiferromagnets 
(AFs) at T > are discussed. The specific case of an impurity coupled symmetrically to two neigh- 
boring host spins is considered. The specific feature of this problem is that the defect is degenerate 
(frustrated) being located in zero molecular field. It is shown that this problem can be described by 
spin-boson model without tunneling term and with a more complex interaction. We demonstrate 
that the effect of the host system on the defect is completely described by the spectral function. It 
is found within the spin-wave approximation that for not too small ui the spectral function is pro- 
portional to uj'^ /J^, where J is the exchange constant between the host spins. The defect dynamical 
susceptibility is derived using Abrikosov's pseudofermion technique and diagrammatic expansion. 
The calculations are performed within the fourth order of the dimensionless coupling parameter /. 
It is found that transverse impurity susceptibility x± {^) has a Lorenz peak with the width propor- 
tional to f^J{T/J)^ which disappears at T = 0, and a non-resonant term. The later term diverges 
logarithmically as uj,T 0. The static susceptibility x(0) has the free-spin- like contribution 1/{4T), 
and a logarithmic correction proportional to ln( J/T). The influence of finite concentration of the 
defects n on the low-temperature properties of AF is also investigated. A logarithmic correction to 
spin- wave velocity of the form nf^ ln\J/uj\ and an anomalous damping of spin waves proportional 
to n/''|a;| are obtained. The results of the present paper can be applied to other systems with a 
frustrated impurity in which the spectral function is proportional to u^. 

PACS numbers: 75.10.Jm, 75.10.Nr, 75.30.Hx, 75.30.Ds 

I. INTRODUCTION 

Defects (or impurities) with intrinsic degrees of freedom exist in many condensed matter systems. Isolated spin in 
metal is the most known example (Kondo impurity). There is a number of other types of such impurities including 
two-level systems in glasses, crystal-field states of the rare-earth ions, degenerate or slightly split Jahn- Teller defects, 
quantum dots etc. The interaction of these impurities with propagating excitations of the host system (electrons in 
metals, phonons, spin-waves etc.) governs the impurity dynamics and the low-temperature thermal and transport 
properties of the host system. 

The widely used model for investigation of two-level defect dynamics is the spin-boson model which Hamiltonian 
has the formic 

Ti-sbm = Hd +^.0 + Hint, (1) 

where the first, the second and the third terms describe, respectively, the isolated defect, the host system and their 
interaction. Then, Hd — — -^Actj; + isCTz, where cr is the Pauli vector describing the defect. Ho is modeled by a set of 
harmonic oscillators: Hq = J2ai^''^a£aQa + ^a/C^'^iQ)] and the interaction term has the form Hint = o'z J2a ^aQa, 
where Cq. are some constants. Essentially, dynamics of the defect is determined by the spectral function characterizing 
the system. Commonly a power-law dependence lu"" of this function is discussed, where u > 0. In the most investigated 
Ohmic case v = 1. Despite its simplicity the spin-boson model has found numerous applications ranging from electron 
transfer to quantum information processing. Meanwhile its modifications are needed in some cases. 

In Refi^ one of us (S.V.M.) has extensively studied the problem of interaction of a defect with intrinsic degrees 
of freedom with 3D acoustic phonons in dielectrics. A specific approach has been proposed in which degeneracy 
of the impurity is assumed to be arbitrary. This approach is based on Abrikosov's pseudofermion technique** and 
diagrammatic expansion. In the case of the two-level defect the Hamiltonian of the model considered in Ref.^ differs 
from Eq. 0J by the absence of Hd (the defect is assumed to be degenerate) and by another type of interaction which 
has the more general form 

H.„t=5E^''^''(^o)' (2) 

A' 

where S ~ ^a, Rq determines the position of the impurity in the crystal, g is the interaction strength, index jj, labels 
Cartesian components, and e^(Ro) are some operators of the host system. It was found in Ref.'^ that similar to the 
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spin-boson model the effect of the host system on the defect is completely described by the spectral function given by 
the imaginary part of the retarded Green's function of operators e^(Ro): 

A^.H = -z/ dte^'^*([e^(Ro,t),e'^(Ro,0)]), (3) 

where (...) denotes the thermal average. In the case of 3D acoustic phonons ImA^i^(tj) is proportional to uj^. The 
approach proposed in Refj^ allows to obtain all the results in a general form independent of the particular view of 
the operators e''(Ro) and the value of S. The only restriction is that the spectral function is proportional to u;^. The 
results for various systems would differ only by some constants. Then, one can discuss the effect of the new terms 
in the interaction in comparison with the spin-boson model. In the case of interaction Q all the components of the 
susceptibility have T-independent non-resonant term and a Lorenz peak with the width F oc where / is the 

dimensionless coupling constant and O is a characteristic energy. The real part of the non-resonant term is a constant 
at |ci;| ^ Q and the imaginary one is proportional to w. At the same time in the spin-boson model the transverse 
susceptibility has only the non-resonant term. It was also shown^ that the scattering on the impurities leads to the 
anomalous 3D acoustic phonon damping proportional to n/'^w^, where n is concentration of the impurities (damping 
caused by scattering on static defects is proportional to w^). Afterward the suggested approach has been successfully 
applied to investigation of defects in glasses^ and in cubic metalsi 

In the present paper we apply the approach discussed in Ref^ with the spectral function proportional to to 
the problem of two-level degenerate defect. As mentioned above, the nature of the defect and the host system is not 
essential. The results will depend on the special form of e^(Ro) in Eq. ||2J) via some constants. Thus, our discussion 
are applicable to all systems with degenerate defect and the spectral function proportional to uj'^. We demonstrate 
below that an example of such a system is 2D Heisenberg antiferromagnet (AF) with the impurity spin-i coupled 
symmetrically to two neighboring host spins (see Fig.^. This is the particular subject of the present investigation. 

Below we show that the spectral function is proportional to if the interaction of the defect with 2D AF is 
determined by spin waves. It is well known that there is no long range order in Heisenberg 2D AF at T > 0^ 
Nevertheless, as it has been shown theoreticallySiiflili and confirmed experimentally-'^^, the spin waves are well defined 
in paramagnetic phase of 2D AF if their wavelength is much smaller than the correlation length ^ oc exp (const /T). 
It is found below that the interaction is determined by spin waves and the spectral function is proportional to 
a u! ^ Ja/£^, where J is the coupling constant between the host spins and a is the lattice constant. Then, a small 
interaction (for definiteness interplane interaction) of the value of ry <C J can stabilize the long range order at finite 
T. It is obtained below that the spectral function is proportional to at w ^ ry for the ordered quasi-2D AF. We 
assume that the interaction of the defect with AF has the form with e''(Ro) 
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Hf„^ =gS(si+S2), (4) 

where Si^2 denote the host spins from different sublattices. For the following consideration the sign of g is insignificant. 

It should be stressed that one must distinguish symmetrically and asymmetrically coupled impurities (see Fig. [T]l. 
Symmetrically coupled impurity is located in the zero molecular field. It remains degenerate and the spectral function 
is proportional to cj^. In the case of asymmetrically coupled impurity, where the molecular field is nonzero, there 
is splitting of the impurity levels and the spectral function has terms with weaker oj-dependence. For instance, we 
demonstrate below that the spectral function for defect coupled to one host spin is proportional to a constant. In this 
paper we consider only the symmetric case. Our results are also valid with certain additional restrictions for slightly 
split nearly symmetrically coupled impurities (see below). 

Previously, different types of impurities in 2D Heisenberg AF have been extensively studied. It is believed that 
this problem has a relevance to the physics of some high-Tj, materials. In such compounds as La2_2;Sr^Cu04 
and YBa2Cu306+2: the mobility of the holes is very small at low level of doping before the onset of supercon- 
ductivity. This finding has generated particular interest to problems of an extra spin coupled to one spin of 2D 
AFj^^i-^'^i^^i-^^i-^'^'i-^^i^^i^^i^-^i^^ and an extra spin coupled to sublattices symmetricallyi22i2i24 It was proposed in Reim^ 
that the appearance of the hole could lead to a ferromagnetic interaction between corresponding two spins of the 
lattice. This work has stimulated studies of a missing or a ferromagnetic bond between one pair of spins in the 
latticeii^iS& As inserting the static non-magnetic defects into the planes by means of replacing of Cu atoms with 
non-magnetic ones (e.g., Zn) is a common method to investigate the properties of Cu02 planes of high- Tc compounds, 
a missing spin in the lattice (vacancy) has been also discussed ii^*2iiS242LSSi2S 

The problem of an added spin in 2D AF at T = has been studied theoretically in Refiii. One of the most 
remarkable findings of that paper is a singular logarithmic frequency behavior of the defect dynamical susceptibility. 
Impurity static magnetic susceptibility x('^ = 0) for 2D AF has been evaluated in Ref>2£ in the case of symmetri- 
cally coupled impurity. It was demonstrated that x(0) has a Curie-like term and a singular logarithmic correction 
proportional to ln( J/T). 
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The defect static magnetic susceptibility in 2D AF being near quantum critical point (QCP) has been discussed 
recently in RefsJi^ using nonlinear sigma model. A classical-like behavior of the form x(0) = S"^ /{2>T) with a 
logarithmic correction proportional to ln(l/T) was obtained. The constants before and under the logarithm were 
found to be universal near QCP, being the same for all types of defects and independent of the strength of the defect 
coupling. It is argued in Rcf.^-'^ that this behavior of static susceptibility also holds far from QCP for vacancy and 
for impurity spin coupled to one host spin when T <C \g\. Meanwhile the constants are no longer universal far from 
QCP. This finding is in agreement with results of numerical simulations^ii^ They have been also confirmed by some 
other theoretical approaches (iS. Both the classical-like form of 1/T-term and the logarithmic correction are related to 
the nontrivial long-range dynamics in the systemiiSiSi 

In the present paper we study Heisenberg 2D AF at T > far from QCP. Our aim is to find the dynamical 
susceptibility of the impurity x{'^) at T > and to discuss the influence of such impurities on the low-T properties 
of 2D AF. These problems have not been addressed yet for symmetrically coupled defects: only the ground state 
properties at T = 0£2i^ and the static susceptibility^ have been discussed. The calculations are performed within 
the order of where / oc g/ J is the dimensionless coupling parameter. We show that the transverse impurity 
susceptibility x^{'^) has a Lorenz peak with the width F cx f'^J{T/J)^ that disappears at T = 0, and a non-resonant 
term. The imaginary part of the non-resonant term is a constant independent of T at 3> and the real part has 
a logarithmic divergence as T ^ 0. Similar logarithmic singularity was found in Ref.^"^ at T = 0. The longitudinal 
susceptibility X||(^) has the non-resonant term which differs from that of x±i^) by a constant and a Lorenz peak. 
We demonstrate that within the order of the width of the peak is zero. Its calculation is out of the scope of this 
paper. 

The static susceptibility has the free-spin-like term S{S + 1)/{3T) and a correction proportional to /^ln(J/r). We 
point out here the sharp difference between symmetrically and asymmetrically coupled impurities that takes place 
in the regime T <C l^l (by asymmetrically coupled impurities we mean here either the added spin coupled to one 
host spin or the vacancy which is the particular case of the added spin with g ^ oo). The leading 1/T-term has the 
free-spin-like form in the symmetric case and the classical-like form in the asymmetric one. Moreover, the logarithmic 
correction is proportional to g^ in the symmetric case and it does not depend on g in the asymmetric one^S*2i The 
difference is related to the fact that the impurity spin coupled asymmetrically aligns with the local Neel order^iSi^i 
whereas the symmetrically coupled impurity is located in the zero molecular field. 

The fact that the spectral function in 2D AF is proportional to only at w ^ {t; or Ja/S,} leads to the following 
restriction on the range of validity of the results obtained: max{T, |a;|} ^ {77 or Ja/£,}. If the defect is slightly split 
(for definiteness by magnetic field H) this condition turns into max{T, ^ max{{77 or Ja/^}, gfisHS}. For nearly 
symmetrically coupled impurity one has: maxj/^, ^ niax{{?7 or Ja/i^}, |gi~52|}, where gi_2 are values of coupling 
with the corresponding sublattices (see Fig. 

The results described above are valid for isotropic interaction We also consider interaction containing only 
one term: Hint = gS^{sf -|- s|). In this case the xx-component of the impurity susceptibility is zero whereas yy- 
and zz- ones have only the non-resonant term. This model is identical to the spin-boson model without Hd- 
The Hamiltonian can be diagonalized exactly and an exact expression for xi^) can be obtained. Below we perform 
the corresponding calculations for the spectral function proportional to and confirm the results obtained within 
our approach. One of the most interesting features of the exact result is that the static susceptibility has the form 
X(0) oc r-i-«, where C oc pT/J. Within the first order of p one has 1 / (4T)-term and the logarithmic correction. 
Thus, we see that in the modified spin-boson model taking into account the higher order logarithmic corrections leads 
to the non-trivial power-law T-dependence of x(0)- 

The influence of the finite concentration n of the defects on the low-temperature properties of 2D AF is also 
considered. For not too small cu we find the logarithmic correction to the spin- wave velocity of the form n/^ln \J/uj\ 
and an anomalous damping of the spin- waves proportional to n/*|a;|. Similar logarithmic correction to the velocity 
and damping were obtained in Ref.^^, where vacancies in 2D AF were studied. It is demonstrated that interaction 
of the spin waves with defects modifies the spectral function which acquires new terms proportional to n exhibiting 
weaker w-dependence. These terms should be taken into account at small enough lu and the problem should be solved 
self-consistently. The corresponding consideration is out of the scope of this paper. Within the range of validity of 
our study we do not obtain a renormalization of the magnetic specific heat which is proportional to in 2D AF 
without impurities. At the same time it was obtained^-'^- that vacancies give rise to a constant contribution to the 
density of states that in turn leads to a large correction to the specific heat proportional to nT. 

The rest of the paper is organized as follows. The model, Abrikosov's pseudofermion and diagrammatic techniques 
employed for the calculations are discussed in Sec.^] The pseudofermion Green's function, the pseudofermion vertex 
and the impurity dynamical susceptibility are derived in Sees. IIII Al UlTO Another type of interaction of the defect 
with the host system are discussed in Sec. IIII Dl The exactly solvable spin-boson model (without Tid) which is a 
special case of our model is also studied in Sec. IIII Dl and a comparison with our results is made. Influence of the 
defects on low-temperature properties of 2D AF is considered in Sec. IIVI The spin- wave spectrum and the specific 
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heat are studied in detail in Sec. IIVI Section Ivl contains our conclusions. A few appendices are included with details 
of the calculations. 

II. MODEL AND TECHNIQUE 
A. Model 

Let us formulate in somewhat detail the model discussed in this paper. We consider systems which Hamiltonian 
can be represented in the following form: 

^ = ( "kOfk + ^ ) + Hint, (5) 



k 



where the first term describes non-interacting low-energy propagating modes of the host system (e.g., phonons or 
magnons) and the second one is the coupling of the degenerate impurity with the system for which we have the 
general expression 0. It is also supposed in this paper that the imaginary part of the function A^p/(w) given by 
Eq. Q has the form 

ImA^^(w) ^ ~A sgn{uj)A{uj)df,^, (6) 

where ^ is a positive constant which dimensionality is inverse energy, Q is the characteristic energy, A(a;) is a cut-off 
function which is equal to unity at \uj\ < and decreases rapidly to zero outside this interval and c?^^ is a tensor. 
As was also mentioned above, the particular nature of the defect and the host system are not essential in our study. 
We will use the general expressions ^ and © in all calculations. Therefore results for different systems would differ 
only by some constants. Nevertheless we discuss now the specific system, Heisenberg 2D AF, which can be described 
by this model. 

It is well known that 2D AF at T 7^ has no long range order.2i The average z-component of the spin in AF is 
given by 



1 4s J - £k 4s J ^ jV(ek) , . 



iV^ 2ek N ^ 

where N is the number of spins in the lattice, s and J are values of the spin and the exchange, respectively, ek is 
the spin-wave energy which is equal to V&sJk at small k and iV(ek) = (e'^''/-^ — 1)~^. The first term in Eq. Q gives 
the well known correction to the average spin at T = which is approximately equal to 0.2. The last term describes 
the spin reduction due to thermal fiuctuations. At ek ^ T we have N{ek) ~ T/ck. Thus the last term diverges 
logarithmically at small k in 2D AF. A weak interaction of the value of 77 ^ J such as anisotropy or an interplane 
interaction can screen this divergence and stabilize the long range order. For definiteness we consider interplane 
interaction. It leads the momentum to become a 3D vector. As a result an additional term appears in the spin-wave 
energy proportional to rik± at small k±_ , where kx is the component of the momentum perpendicular to the plane of 
the lattice. As a result the last term in Eq. |7J is small and the spin waves are well defined if 

-\n( -)<$:!. (8) 



S J \ ST] 

We will assume below that this condition holds. 

It is shown in Appendix 1X1 that within the spin- wave approximation the function ImA^i,(cj) has the form © for 2D 
Heisenberg AF when |a;| ^ {77 or Ja/^}. The particular expressions for Q, A, d^v are also established in Appendix 1X1 
Within the spin-wave approximation the only nonzero components of are xx- and yy- ones provided that z-axis 
is directed along magnetization of the sublattices. 

Abrikosov's pseudofermion representation of the impurity spin S is used below. The value of S is assumed to be 
arbitrary in this approach. Nevertheless we restrict ourself in this paper by S* = 1/2. It is demonstrated below that 
the matrix structure of the pseudofermion Green's function and the vertex is much simpler in this case. Consideration 
of larger impurity spins is out of the scope of the present paper. 

We point out that a new approach has been suggested recently in Refsi^i*^ for spin-i impurity problem. This 
approach bases on Major ana- fermion representation of the impurity spin. It was demonstrated that this representation 
simplifies significantly analysis of the model if one can restrict ourselves by first terms in the expansion by coupling 
parameter. It will be clear soon that in our case the question of possibility of such restriction requires analysis of the 
diagrams of the third order within this approach. Carrying out of such analysis is out of the scope of the present 
paper. 
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B. Abrikosov's pseudofermion technique 



We use below Abrikosov's pseudofermion technique for the calculation of the impurity dynamical susceptibility. It 
was suggested in Refi^ for Kondo effect investigation (see also Refsi^2i24 for discussions). The same approach has 
been applied for the problem of impurity in other systems by one of us (S.V.M.) in Refsi2*Sii. Let us formulate this 
technique briefly in the convenient for our purpose form. 

The impurity spin S is represented as 



(9) 



where m is the spin projections, a]^j and Um are operator of creation and annihilation of some particles (fermions for 
definiteness). It is easy to verify that the spin commutation rules are satisfied in this representation. A wave function 
of the impurity is characterized now by the occupation numbers of 25" + 1 states: \ns, ns-i, ■ ■ ■ n-s)- Obviously, the 
states with zero or more than one particles are not physical ones and we have to eliminate them carrying out the 
thermodynamic average. As a result the thermodynamic average of some operator Y has the form: 



Y 



TrP'^y'ipY) 



(10) 



where the traces are limited to physical states and p — exp(— 7i/T) is the statistical operator with the Hamiltonian 
Ti. given by Eq. (Q. This representation is not convenient not allowing to use the standard diagrammatic technique. 
To overcome this obstacle an additional term in the Hamiltonian is addedi^ 



XN, 



pf 



(11) 



where Npf is the number of pseudofermions. We show now that the average of Y has the following form equivalent 
to lfTn|: 



Y = lim 



_ Tr(py) 

A^oo Tl^pNpf) ' 



(12) 



where p = exp{— 7Y/T} and H = H + Tix- We consider in this paper such Y that do not contain terms without 
pseudofermion operators. As a result there are no contributions both to numerator and denominator in Eq. H12|l from 
states with no particles. As Ti. does not change the number of particles in state, we have for the matrix elements: 
Pki = Pki exp(— A'^;A/T), where Ni is the number of pseudofermions in states \k) and Therefore, contributions from 
the states with more than one particles are exponentially small compared to those from the physical states which are 
proportional to e~^l'^ . The common factors e~^^^ in the numerator and the denominator of Eq. (|12|l cancel each 
other. Then contributions from states with one pseudofermion survive only in the limit of A ^ cx) and Eqs. H1U|I and 
(|12|l appear to be equivalent. 

Quantities in the right part of Eq. H12|l can be calculated using the conventional diagrammatic technique. The 
Hamiltonian in the pseudofermion representation has the formi^ 



^ = (E [^^k + + aj„am ) + .9 E «L'^rn'™«rne^(Ro) = + 



(13) 



The dynamical susceptibility of the impurity in the representation of interaction can be written using Eqs. H12I) and 
(|13|l in the following formi^ 



XpiiuJn) = lim TV 



A — >oo 



1/T 



dTe*""^Tr 



e~^°/^Tr (p(r)P(0)e ( l 



T 



'HIT 



E"^ 



(14) 
(15) 

P = a'^Pm.m'a-m', (16) 

mm' 

where P is a spin projection. In the zeroth order of the interaction l-Lint we have M — {2S + l)e^^^^ . 
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C. Diagrammatic technique 

First diagrams for xp{^) f^nd a graphical representation of the result of all diagrams summation are shown in Fig. 12 
where thin lines with arrows represent the bare particle Green's functions: 

= (17) 

and wavy lines denote bosons Green's functions A^^/(ia;„). Notice that the diagrams with only one pseudofermion 
loop should be taken into account because, as is seen from Eq. (|14|l . each loop is proportional to the small factor of 
e~^l^ . The contribution from diagrams with one loop is finite because their factor of e~^l'^ is canceled by that from 
M. 

For the calculation of Xp('-^) "^6 use below the diagrammatic technique employed in Refsi^iSii. Let us discuss it 
briefly. Firstly, we have to make an analytical continuation of diagrams for the dressed pseudofermion Green's function 
and for the vertex Vp(lo\ + w, wi) from imaginary frequencies to real ones. Then, we have to express xp(a;) via these 
quantities. 

To make the first step of this program let us choose frequencies of wavy lines to be independent variables over 
which the summations are taken. It can be done in such a way that these frequencies are contained in arguments of 
G'"^-functions with positive sign (see Fig.|2J). Then, each sum over a discrete frequency can be replaced by an integral 
over a contour enveloped the imaginary axis with an additional factor of (27ri)~^A^(aj): 

where N{uj) — (e"/"^ — is the Plank functioni^iSiLS^ The contour can be deformed so as to embrace the real axis. 
In evaluation of the resultant integral one should not take into account poles of G^^-' -functions because residues in 
these poles are proportional to N{X) « exp(— A/T). At the same time functions A^p/(a;) has a discontinuity on the 
real axis equal to 2iImA^^/(a;). As a result all contour integrals can be easily transformed to those over the real axis 
and we lead to the following diagrammatic technique: each wavy line corresponds to 7r~'^iV(u;)ImA^^' (w); frequencies 
of wavy lines should be taken so as they are contained in arguments of G^*^-* -functions with positive sign; integration 
over all frequencies of wavy lines is taken in the interval (—00,00). 

One can conclude from analysis of all concrete diagrams for the vertex Vp{iu}i + iu),iuji) that it is an analytical 
function of two independent variables iuji + ioj and iwi with cuts along the real axis. A general proof of this statement 
has been also given«25i 

As a result we have for the dynamical susceptibility after the analytical continuation from the discrete frequencies 
to the real axisi^iSiL2& 

/oo 
da;e-^/^Tr{F [G{x + u)T++{x + w, x)G{x) - G*{x)rp- {x, x - lu)G*{x - w) 
-00 

- G{x + Lu)T+- {x + UJ, x)G* (x) + G{x)r+-~ {x, X - Lo)G* {x - Lo)]}, (19) 

where G{iu) is the retarded Green's function, the trace is taken over projections of the impurity spin and signs at 
superscript of Fp denote those of imaginary parts of the corresponding arguments (e.g., Fp~(a;, y) = Tp{x + iS,y — 
i6)). An energy shift by A has been performed during the derivation of Eq. (|19|l . As a result the Fermi function 
(g(a:-i-A)/T _|_ lias been replaced by exp(— (a; -I- A)/r) and the functions G and Fp no longer depend on A. These are 
those functions we calculate in the next section by the diagrammatic technique. It is clear that the bare pseudofermion 
Green's functions in this case are G|^^/(a;) — 5mm' 1^- 



III. DYNAMICAL SUSCEPTIBILITY OF THE IMPURITY 



We derive analytical expressions for the dynamical susceptibility of the impurity in this section. Perturbation theory 
is used for this purpose. It can be done if the dimensionless constant 

= ^ (20) 

is small. Meanwhile we have to take into account also terms of the order of because the finite width of the Lorenz 
peak in the dynamical susceptibility arises in this order. In 2D AF we have from Eqs. (|A10|) : p = {s/n / As){g / JY . 
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A. Pseudofermion Green's function 



We turn to the calculation of the pseudofermion Green's function Gmm'i'jj)- The Dyson equation for it has the 
following form: 



^Gmm'{^) — ^mm' + / ^ 2]mn(w)G„m' (w). (21) 



The first diagrams for Sm„(cj) are presented in Fig.Ol Let us discuss its matrix structure first. It is determined by 
corresponding products of operators S'^ and tensors d^^- For example, this combination for the second digram in 
Fig. 13 has the form S^S'^S'^ S'^ d^,fi'dyui- As is shown in Appendix IbI all such combinations are proportional to the 
unit matrix for the two-level impurity. It should be pointed out that there is no such simplification in the case of the 
impurity with the value of spin greater than 1/2. As a result the equations for the Green's function and the vertex 
become more complicated. The corresponding consideration of the large-spin impurities is out of the scope of the 
present paper. 

Taking into account its matrix structure we have for the Green's function: 

G„im'{^) = 5„„n>G{uj) = (22) 

ijj — Zj(u;) 

The diagram of the first order shown in Fig. Ogives the following contribution to '^{lo): 

= r dx^^^N{x)K{x), (23) 

J-oo X + OJ + 10 

i?s = S^'S'^d,,, (24) 

where the constant p is given by Eq. H2Q(I and i?s = 1/2 for 2D AF. It is convenient to extract from this expression 
a term proportional to uj as follows: 

I](i'(w) = / dx LJ N{x)A{x)-Rs^ dxxA{x), (25) 

J-oo X + Lu + id 7r6 Jq 

where we have used that A{uj) is an even function and N{—x) = — 1 — N{x). Notice that the second term in Eq. H25() 
is the T-independent constant. Then it can be included in the renormalization of A and omitted. The first term 
in Eq. H25|l is proportional to f'^ujT\n{T/u!) at small w. It would seem that a great renormalization of the Green's 
function takes place if f'^T\n\T/uj\ > Q. Meanwhile we show now that the logarithmic singularities at real uj are 
screened by a finite damping which is of the order of /"*. Let us represent the Green's function in the form 

Gi.) = 1-^, (26) 

where Z{uj) and ^{uj) are some functions, 7(1^') is a real one, and a constant term in the denominator has been 
attributed to the renormalization of A and discarded. Evaluating contribution from the first diagram in Fig. |3| using 
Eq. 1)26(1 we have in the first order a correction to the constant. 



Z^'Huj) = R^^ / dx 



X + UJ + i"f{x + uj) 



N{x)A{x) (27) 



and 7(^)(w) = 0. The logarithmic divergence in expression H27fl at real uj is screened by the term ij{uj + a;) in the 
denominator. It is shown below that 7(w) is proportional to f^T^ at |w| ^ T. Contributions to Z{uj) from the higher 
order diagrams are also small by the same reason and we can restrict ourself by the first correction to it. 

To obtain 7(0;) one has to take into account /^-terms from the first and the second diagrams shown in Fig.O Along 
with the small corrections to the constant and to Z{uj) we have: 



7(w) = -Rj^T^ I dx\x{x+uj)\N{x)N{-x~ui)A{x)A{x + uj), (28) 



— 00 



R^ = S^S^iS" ,S^']df,^,d,,,. (29) 
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It is important to note that 7(w) is the constant at \uj\ <^ T: 



^efl) ^ ,30, 



e 



Notice that from the physical reason and /q should be positive. For instance, R^ = 1/4 for 2D AF. 

It is significant to note that ImZ(cj) and 7(0;) have the following property at 3> /q which will be useful in the 
following: 



7(-w) = e-"/^7(w). 
In fact these functions are exponentially small at negative w if 3> T". 



ImZ(-cj) = -e--/^ImZ(a>), ^^^^ 



B. Pseudofermion vertex 



Let us turn to the consideration of the pseudofermion vertex Tp{x + w,a;). First diagrams for this quantity are 
presented in Fig. 01 It is shown in Appendix IbI that T Pmm' (x + w, x) is proportional to Pmm' ■ Thus, it is convenient 
to introduce a new quantity: 



J., , ^ PTp{x + uj,x) 

r{x + uj,x)^ = , (32) 

where we use the following notification: Y — Tr(Y). As is seen from Eq. (|19|l . we need four different branches of 
T{x + uj,x). It is clear that r++ = (F^^)* and within the first order of one has: 

r++{x + u;,x) ^ 1 + Ri^ dyy\y\N{y)Aiy)Gix + y + cj)Gix + y), (33) 



R, ^ (34) 

It is seen from Eq. (|33|l that the poles of G-functions in the integrand are on the one hand from the real axis. Hence, 
the second term in Eq. H33|l is much smaller than unity and we can restrict ourself by this precision. 

The situation is different in the case of = (F ^)*. The first correction to it is given by Eq. H33|l with G*{x + y) 
put instead of G{x + y). Therefore, poles of the Green's functions are on the opposite sides of the real axis. As a 
result at a; = the integral diverges at finite x as /q ^ and one has to sum all series to determine F^ . We write 
now an equation for F^ in which the most singular diagrams in each order of are taken into account. As a result 
of analysis of the diagrams up to the fourth order of we have obtained that in the most singular diagrams each 
wavy line connects points from different sides of the vertex (like in the second, the third and the fourth diagrams in 
Fig.^land not like in the last one) and crosses no more than one another wavy line. Thus, the most singular diagrams 
are taken into account in the following equation: 

F+~(a; + w,x) = 1 + Ri^ / dyy\y\N {y)K{y)T+~ {x + y + lo,x + y)G{x + y + lo)G* {x + y) 



TT 



e 



—00 
2 



+■^2 ( ^ ) / rfyic?2/2yiy2|yiy2| A^(l/i)A^(2/2)A(yi)A(y2)r+ {x + yi+y2 + uj,x + yi+ 2/2) 



Tre^ 

X G{x + yi + uj)G{x + yi+y2+ uj)G*{x + yi+ y2)G*{x + y2), (35) 



R2 = = ^ . (36) 

The second and the third terms in Eq. H35|l take into account diagrams with a rung and with crossing of two neighboring 
rungs, respectively. Let us try to solve this equation by iterations. It is easy to verify that at \uj\ ~ Pq and Jq — » the 
divergence in the second term occurs in the second iteration only. Moreover it is of the same order as the divergence 
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of the third term in the first iteration. As a result the equation H35|l can be rewritten as follows: 

j-2 roc 

r+-ix + u;,x) = l + / dyy\y\N{y)Aiy)G{x + y + iu)G*{x + y) 

+ (:^) / dyidy2yiy2\yiy2\N{yi)N{y2)Aiyi)Aiy2)T^^ (x + yi +y2 + uj,x + yi +^2) 
X [RlG*{x + yi)+R2G*{x + y2)]G{x + yi+uj)G(x + yi+y2+uj)G*{x + yi+y2). (37) 



It can be solved if one notes that at \uj\ ^ T, 8 in the third term the area of integration near poles of two last Green's 
functions is essential. As the rest factors of the integrand changes slightly at such j/2 one can set in them y2 = —x — yi 
and neglect their dependence on lu. As a result we obtain the following equation: 

r+-{x + Lj,x) = l + / dyy\y\N{y)A{y)G{x + y + u;)G*{x + y) + K{x)r+-{cj,0)—-—, (38) 

- {01 ^^ t'^'^t!;::';, - ^)-^"')^'- ^')- <-> 

It can be easily solved with the result: 

r+-{x + u;,x) ^ 1 + dyy\y\N{y)A{y)G{x + y + u;)G*{x + y) + K{x)Zr——-, (40) 

ttB uj + 2ir 

r = ro-7rK{0)^27rRri^] dyy' N iy)[l + N {y)]A^ (y) , (41) 

j-2 /•CO 

Zr = l + i?i^/ dyy|y|iV(y)A(y)|G(y)p, (42) 



Rr = R,-RUR2=''''''^^'"^p^''^'-''^r (43) 

In 2D AF Rr = 1/2 and for P = S'^, S'^ and P = S^', respectively. Evidently, the temperature dependence of F for 
P = , S'^ is the same as that of Pq given by Eq. H3U|I. Note, the third term in Eq. H4U|I is much smaller than the 
second one when ::§> Pq/P- 



C. Properties of the impurity susceptibility 

We can derive now the impurity susceptibility using the general expression (|19|l . Eqs. (|26|l . (|27|l . (|28|l and H30(l 
for the Green's function and Eqs. 133|l and 1)40(1 for the branches of the vertex. As a result of tedious but simple 
calculations presented in Appendix El we have for the dynamical susceptibility of the impurity up to terms of the 
order of f^- 



\x\N{x)A{x) 
x + 2iPo 

2^y^ / dx (. (44) 



{ 2^^ / ^r. r , \x\N{x)A{x)\ 2iPo ^ p , 

2 [ P{uj + 2iP) \ ttQ X + 2iPq J T(w + 2iPo) ttO 

P r , sgn{x)A{x) 



i^. ^i^.-i^.^ ^^^y]-^- . (45) 

The first term in Eq. (|44|l is the Lorenz peak with the width P. The second one is a Lorenz peak with the width 
Iq and small, proportional to p, amplitude. The last term is the non- resonant part of the susceptibility. The 
imaginary part of the non-resonant term in Eq. ((44|l at Iwl S> Pq is proportional to sgn(ci;) and the real one contains 
the logarithmic singularity of the form ln(Li;^ + Pq). At T — and uj the nonresonant contribution survives only 
and the susceptibility has the logarithmic singularity. Such a singularity has been obtained for the two-level impurity 
at T = in Refiii. The first and the second terms in Eq. H44|l are calculated assuming that <C T. At ^ T 
these terms are of the order of and their taking into account exceeds the range of accuracy. 
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It should be noted once more that the particular nature of the defect and the host system is not essential in the 
above consideration. We use the general expression ||2Jl for the interaction and assume that the function ImA^i^((jj) 
has the form . Apart from only coefficients R in the resultant expression 1)44(1 depend on the nature of the defect 
and the host system. We demonstrate now that systems with different symmetries of the interaction show different 
behavior of the impurity. In the case of isotropic interaction all the components of tensor d in Eq. (jSJ are nonzero and 
all components of the impurity susceptibility have the same structure: the Lorenz peaks and the nonresonant term. 
When one of the component of d is zero, say zz- one like in 2D AF, the behavior of transverse components Xxi^) 
and Xy{^) differs from that of longitudinal one Xz(w). The transverse components have the Lorenz peak -T ~ and 
the nonresonant term. The Lorenz peak with the width Iq disappears because Ri — 0. The transverse component 
contains the nonresonant term but F = (see Eqs. (I41II and H43|l ') and our precision is insufficient to determine the 
resonance terms in Xzi'^)- The corresponding calculations of Fq and F with higher precision is out of the scope of 
the present paper. If only dxx is nonzero, then Xxii^) = whereas Xyi'-^) Xzi'-^) have only the nonresonant term. 
This particular situation is considered in detail in the next section. 

For static susceptibility xp(0) we have from Eq. (|^ : 

XP(0) - + (46) 



WiT) = 2R^^ r dx^i^^MM (T - xiV(x)) , (47) 



where the first and the second terms in W{F) stem from the non-resonant and the resonant parts in Eq. 144|l . 
respectively. It is easy to verify that W{F) is a real value up to the order of (remember, the susceptibility has been 
calculated up to terms of this order). Integrations in Eq. l|T7jl can be simply carried out and we have at T <C 6: 



Xp(0) = 5(l-2i?.^)+P^2i?,^ln(|). (48) 



Thus the uniform susceptibility has the free-spin-like term P^{2F)^^ which amplitude is slightly reduced by the 
interaction and the correction proportional to /^ln(8/r). Expression (|48|l is in accordance with that of RefiSi. 
Similar result for the static susceptibility, 1/T-term and a logarithmic singular correction to it, has been obtained in 
2D AF near QCPi2i It was found that the static susceptibility exhibits the classical-like Curie behavior of the form 
S"^ / {2tT) and the coefficients before and under the logarithm are universal values independent of the particular type 
of the impurity and the strength of its coupling to the host system. Remarkably, this behavior remains also far from 
QCP for asymmetrically coupled impurities (vacancy and added spin) at T <C g| i3^i4^i6,2i ^j^g constant under 
the logarithm becomes non-universal. These findings are related to the nontrivial long-range dynamics of the 2D 
AF. Then we point out the significant difference between dynamical properties of symmetrically and asymmetrically 
coupled impurities in the regime T <^ g: the leading 1/T-terms have the free-spin-like behavior and the classical-like 
one, respectively. Moreover, the logarithmic corrections is proportional to in the symmetric case and it does not 
depend on g in asymmetric one^SiSi As was also pointed out in Introduction, the difference can be explained by 
the fact that the impurity spin coupled asymmetrically aligns with the local Neel orderiiSiSi aX F <^ g whereas the 
symmetrically coupled impurity is located in the zero molecular field. 

It is convenient from this point on to neglect in Eq. 144|) the small corrections and use the following simple expression 
for the transverse susceptibility x^{'^) = Xx{'^) = Xy{^)'- 



Fiuj + 2iF) ttB x + lo + iFo 

iT^xd^) = ^ V(u;2 + 4j.2) +P^Rx-Q^{^)^Z<^)- (50) 



The imaginary part of the non-resonant term is calculated at |cj| Fq in Eq. (|50() . We see that it does not depend 
on the temperature at such lo. It is seen from Eq. (|49|) that the non-resonant term gives the main contribution to the 
susceptibility when 

|a;| » C.0 = /'e ("I") . (51) 



As the sign of Imxp(a;) should coincide with that of lo, the constant R^ given by Eq. 145|) should be positive. For 
example, we have for 2D AF: R^ = 1/2. 
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It should be stressed that there is a restriction on the range of vahdity of the resultant expressions H44() and l|49() 
for Xp(^) ill the case of 2D AF. It is the consequence of the fact that the function ImA{a;) has the form © if 
uj ^ {f] oi Ja/S^} only. It is easy to see that in all calculations performed above one can use the function of the form 
© if the following condition on uj and Fq holds: maxjlo, ^ {77 or Ja/£^}. 



D. Comparison with the exactly solvable model 

Let us consider the special case when the interaction contains only one term: Timt — gS^ e'^CRo). It is seen from 
Eqs. 129() and (|43|l that F ^ Fq ~ and the resonant terms are zero at w 7^ in Eq. (|44|l . One can conclude from 
Eqs. (jUJ and 1(13 that Xa:(w) = and 



P 



(52) 



As it was pointed out above, our model in this case is equivalent to the spin-boson model ^ without Ti^. The 
corresponding Hamiltonian can be diagonalized exactly and an exact expression for the impurity susceptibility can 
be derived. We perform in this subsection the corresponding calculations of x'-^-'(cij) and confirm our results obtained 
above. 

The detailed discussion of the exactly solvable spin-boson model is also necessary by the following reason. It was 
important for our consideration that in expression H26(l for the Green's function is nonzero, i.e., the coefficient 
R~f given by Eq. (|29|l is finite. As is demonstrated in Sec. IIII K\ the imaginary term ij{co) in the denominator of the 
Green's function G{uj) screens the logarithmic singularity of the self-energy part allowing us to restrict ourself by the 
second order of p. In the opposite case, when imaginary part in the denominator of G{lu) is zero, one has to sum all 
the series to calculate the Green's function in the region of lu determined by the condition 



> 1. 



(53) 



It should be stressed that Eq. H52|l for x'^^H'^) valid if the condition (|53|l on the frequency uj is not fulfilled. If the 
condition (|53ll is fulfilled the results of the exact diagonalization should be discussed. 

We represent Hamiltonian of the spin-boson model describing degenerate defect in the following form: 



n^l E(^kQkO-k + PkP-k) + gs^ 



k 



(54) 



where the symbols e^, Qk and Pk stand for the frequency, normal coordinate and momentum of the system propagating 
modes (bosons) with momentum k, where [Qk, P-k'] = "iSkw , Qu = ak + aLkJ -^k = — *ek(Q^k — Q^Lk)- The last term 
in Eq. H54|) describes coupling of the impurity with the system, where -fk is a coupling parameter. In this problem 
definition the spectral function given by Eq. (PJ has the form 



ImA(w) 



2^ 

k 



l^k^ 
Ck 



[S{u - Ek) - S{u! + ek)]. 



(55) 



We treat the elementary excitations within Debye approximation and assume that the spectrum is linear in k: Ck = ck. 
We make also one more assumption which do not effect the results: k is a two-dimensional vector. In this case the 
spectral function is proportional to when the coupling parameter is linear in k: 



1 



Ck, 



(56) 



where V is the volume of the crystal and 8 is the cut-off frequency In this case the spectral function H55|) has the 
form © with A = (4c^)"^ 

The Hamiltonian 154(1 can be diagonalized exactly. It is convenient for this purpose to apply the following canonical 
transformation:'^'' 



k k 



R 



(57) 
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It can be shown that the correlation function can be brought to the form: 



.30 



l] + [iV(6k) + l][e- 



-iekt 



!]}• 



As a result we have for the transverse dynamical susceptibility using the representation I{t) = X{t) — iU{t): 

1 



One obtains for X{t) and U{t) from Eq. 

X{t) = 
U{t) = 



P 



dx[cos(a;i) - 1][1 + 2N{x)], 



dx sin(a:i) = 



P 1 ~ cos(et) 

TT et 



(58) 
(59) 

(60) 

(61) 
(62) 



where p is given by Eq. (jJOJ. It is seen from Eq. H62(l that U{t) is a bounded function and we have in Eq. H60|l : 
sin[J7(t)] w [^(0- The function X{t) along with negligibly small terms of the order of p has another one which is 
large at tT > exp{e/(/2r)} and for which we have with the logarithmic precision: X{t) « -2pT{-KQ)~^ \n{tT)^ 
As a result we obtain: 



P 



2tt \T 



1 - cos(ie) 



(^e)i+c 







'1 


H 


LU 


c 


1 


1 - 


UJ 


c 




C" 






2 


" e 


H 


^2 


e 




e 









'1 


U 




c 


1 


1 - 


UJ 






2 


" e 




2 


6 



sgn(8 — a;) — 



LU 


c" 


e 





sgn(Lj), 
(63) 



where C = 2/2T(7re)-^ We recover Eq. lO from Eq. ^ if the condition ((23 is not fulfilled and \lu\ < 8 or 
8. At w ^ 8 Eq. H63|l can be represented in the simple form: 



C r 



(64) 



At small enough uj, when H53|l holds, equation H52I) is incorrect and we see from Eq. H64|l that the susceptibility shows 
the nontrivial cu- and T- dependences. In particular, the static susceptibility is proportional to T^^^'». In the order 
of p the static susceptibility has the same structure as that obtained above (see Eq. (^HJ): it has the conventional 
term 1/(4T) and the logarithmic correction to this term. Thus, we see that in the modified spin-boson model taking 
into account higher order logarithmic corrections results in the non-trivial power-law T-dependence of x(0). 

IV. INFLUENCE OF THE DEFECTS ON THE HOST SYSTEM 

We discuss in this section the influence of the defects with finite concentration n on the low-temperature properties 
of 2D AF. The spin-wave spectrum and the specific heat of AF are considered below in detail. It will be assumed 
that n ^ 1 in order to neglect interaction between impurities. 

A. Spin-wave spectrum 

The spin-wave spectrum is determined by the poles of the spin Green's functions. The Green's functions of 2D AF 
with impurities are investigated in Appendix^ It is demonstrated there that their denominator has the form 



V{uj, k)^uj^-el + As^g^nx±{iu) [ Jo - Jk cos(kRi2)] , 



(65) 
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where R12 is the vector connected host spins coupled to impurity (it is assumed for beginning that this vector is the 
same for all defects), Jk — 2J(cosfc^ + cosfc^), where J is the exchange constant, and ek = sy^ Jq ~ J^- 

Let us discuss the spin- wave spectrum near points k = and k = ko, where ko is the antiferromagnetic vector. It is 
seen that expression ()65f) is symmetric under replacement of k by k ± kg. Thus, we consider below only the vicinity 
of the point k = 0. We have from Eq. (|65|l : 



2 



l-^e4k)x^(a;) 

ZTT 



(66) 



1 (kR 



(k) = ;^ + ^^^, (67) 



^2 



where it is used that the unperturbed spectrum is linear at A: <C feg: ek = cfc = ^/SsJk. It is seen from Ec^s. I|66() 
and (|67|l that the spectrum appears to be dependent on the direction of the momentum k as a result of interaction 
of magnons with the defects. This circumstance is a consequence of our assumption that the vector R12 is the same 
for all impurities. In fact, it can have four directions and the value (Ri2k)^/A:^ can have two different values: cos^ (/)k 
and sin'^ (/)k, where 0k is the azimuthal angle of k. It easy to realize that M(k) = 1 if all four ways of coupling of the 
impurity with AF are equally possible. 

It is convenient for the following to consider separately the cases of \uj\ ^ loq and \u)\ <C wq, where wq is given 
by Eq. (|51|l . In these cases the non-resonant and the resonant parts, respectively, are dominant in the impurity 
susceptibility (j^ . 

|a;| ^ wq. One obtains from Eqs. (|49|l . (|5()|l and H66|l for the magnon damping 7k and the renormalized spin- wave 
velocity Ck: 

7k = l^l^«(k), (68) 
^u(k)) , (69) 



;;2 _ „2 



cM 1 - In 



where we take into account that i?^ — 1/2. It is seen that the interaction with the defects leads to strong damping 
which is proportional to uj and to the logarithmic correction to the spin-wave velocity. It would seem that at small 
enough k the spin- wave velocity becomes imaginary signifying a phase transition in the system. Meanwhile our theory 
is not applicable at such small k. The interaction with the defects changes the function ImA^;^(a;) as well and this 
renormalization is strong at small enough u. Indeed, one has to use renormalized spin Green's function derived in 
Appendix IdI to evaluate ImA^i^(ix'). As a result of simple calculations similar to those presented in Appendix 1X1 we 
obtain that at \ui\ ^ 8 in addition to the term proportional to w^sgn(w) there is another one proportional to sgn(ci;): 

— j sgn(w)d^t. - sgn(cj) — - — d^^, (70) 
B^-!^[ ,,(l+cos(kR,2))-(Jo-./k)- ^ 

where the integral is taken over the chemical Brillouin zone. The first term in Eq. (|70|l is greater than the second one 
when 

|w| > 0.0267^/2. (72) 

This condition determines the range of validity of our theory at |a;| wo- At w given by (|72|l the logarithmic 
correction to the spin-wave velocity in Eq. (|69|l is small. 

\uj\ <§; Wo- In this case the impurity susceptibility (|49(l is determined by the resonant term. We have for the spin- wave 
damping and the spin- wave velocity from Eqs. H49|) . (|50|l and (|66f) : 



T uj^ + An 27r 

As in the case of 3> ti^o, one has to take into account the renormalization of the function ImA^i/(aj). After simple 
calculations we obtain at \uo\ ^ <d: 

ImA^^(w) = -A sgn{uj)df,^ - 25n/^B ^^^^^^^^^^ d^,^, (75) 
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where the constant B is given by Eq. H71|l . As a result the range of vaHdity of our consideration is determined by 

M<^^»0^004„/^^^ (76, 

It is easy to show that the spin-wave damping H73|) and the correction to the spin-wave velocity in Eq. 174|) is small 
at such to. 

It should be noted once more that if one of the conditions (|72() or H76|) is violated the interaction between defects 
becomes important and our approach is wrong. 

B. Specific heat 

We proceed with the discussion of the magnetic part of the specific heat C (T) . It is convenient to use the following 
formula for its evaluation: 

where the first term describes magnetic excitations and index i in the second term labels impurities. To calculate the 
first term in Eq. (|77|l we use bilinear part of the Hamiltonian (jA5|l and spin Green's functions (|D2|) renormalized by 
interaction with impurities. As a result of simple calculations one obtains up to inessential terms not depending on 
the temperature: 

(4"k + ^ = 51 ^kA^M + 2Arn/2eX y" dLuN{Lj)lmxx{uj), (78) 
X ^'^[ ^i^ -^o-^^kCosCkR,,) ^ ^^^^ ^^g^ 



27r4 



where the integration in Eq. (|79|l is over the chemical Brillouin zone. The diagram for the second term in Eq. (|77|l is 
shown in Fig.|S| It contains the impurity susceptibility calculated above and A^^{lu): 




dujN{uj)[lmx^ioj)ReA^^{uj) + ReXi.(c^)IniAi,^(w)]. (80) 

Evaluating expression H8U|) and summing it with Eq. (|78|l one leads to the following expression: 
F 1 2nf^ f°° f°° 

1^ ^ wE^kA^(^k)-^ / du;N{Lj)u;\u;\A{u;)Rex±{u;)+nfeX du;N{Lj)lmxA^), (81) 

where the constant X is given by Eq. (|79|l . The first term in Eq. H81|) describes the energy of 2D AF without 
impurities. It is equal to Eq cx NT^ /Q^ at T <C 6 and we come to the well known result: C(T) cx for pure 2D 
AF. Evaluation of corrections to Eq from the second and the third terms in Eq. H81(l is a tedious but straightforward 
work. Unfortunately Eq. H49|l for x±{^) has the limited range of validity which is discussed in the previous section. 
Therefore one can not carry out the integration in the second and the third terms of Eq. 181|) in the whole range of uj. 
We have evaluated these terms after the integration over lu at which Eq. 149|l is valid. There are corrections having 
weaker T-dependence than Eq stemming from both the resonant and the nonresonant parts of x±{^)- Meanwhile 
restrictions (|72|) and H7t)|) result in the considered corrections to be bounded below on T and to be smaller than Eq. 
Thus, we do not obtain a renormalization of the specific heat within our precision. 



V. CONCLUSION 



We discuss the dynamical properties of the impurity spin-i in 2D and quasi-2D Heisenbcrg antiferromagnets (AFs) 
at T > 0. The specific case of the impurity that is coupled symmetrically to two neighboring host spins is considered. 
It is shown that this problem is a generalization of the spin-boson model without the tunneling term and with a 
more complex interaction. It is demonstrated that the effect of the host system on the defect is completely described 
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by the spectral function S{uj) which is proportional to oj^. It is found within the spin- wave approximation that the 
spectral function has this w-dependence in 2D AF for not too small uj. In isotropic 2D AF at T > J{oj) oc 
when u) ^ Ja/£,, where J is the exchange constant between the host spins, ^ is the correlation length and a is the 
lattice constant. For ordered 2D AF J^{llj) oc o;^ at o; ^ 77, where 77 <C J is the value of interaction (for definiteness 
interplane interaction) stabilizing the long range order at finite T. 

We stress that one must distinguish symmetrically and asymmetrically coupled impurities (see Fig. ^ . Symmet- 
rically coupled impurity is located in the zero molecular field. It remains degenerate and the spectral function is 
proportional to lu'^. In the case of asymmetrically coupled impurity, where the molecular field is nonzero, there is 
splitting of the impurity levels and the spectral function has terms with weaker w-dependence. For instance, we 
demonstrate that the spectral function for impurity coupled to one host spin is proportional to a constant. In this 
paper we consider only the symmetric case. Our results are also valid with certain additional restrictions for slightly 
split nearly symmetrically coupled impurities (see below). 

The defect dynamical susceptibility x('-^) is derived using Abrikosov's pseudofermion technique and diagrammatic 
expansion. The calculations are performed within the order of where / oc g/J is the dimensionless coupling 
parameter. For our study the sign of / is insignificant. We show that the transverse impurity susceptibility XJ_('-^) has 
a Lorenz peak with the widths F oc f'^J{T/J)^ that disappears at T = 0, and a non-resonant term. The longitudinal 
susceptibility X||('^) the non- resonant term which differs from that of x±{^) by a constant and a Lorenz peak. 
The width of the peak is zero within the order of Its calculation is out of the scope of this paper. The imaginary 
part of non-resonant term is a constant independent of T at |a;| 3> -T and the real part has a logarithmic divergence 
as w, T — > 0. Similar logarithmic singularity was found in Refii at T = 0. 

The static susceptibility has the free-spin-like term S{S + 1)/(3T) and a correction proportional to /^ln(J/r). We 
point out here the sharp difference between symmetrically and asymmetrically coupled impurities that takes place 
in the regime T <^ \g\ (by asymmetrically coupled impurities we mean here either the added spin coupled to one 
host spin or the vacancy which is the particular case of the added spin with g ^ 00). The leading 1/T-term has the 
free-spin-like form in the symmetric case and the classical-like form in the asymmetric one. Moreover, the logarithmic 
correction is proportional to in the symmetric case and it does not depend on g in the asymmetric one^S*2i The 
difference is related to the fact that the impurity spin coupled asymmetrically aligns with the local Neel order^iSiSi 
whereas the symmetrically coupled impurity is located in the zero molecular field. 

The fact that the spectral function in 2D AF is proportional to only at w ^ {77 or Ja/^} leads to the following 
restriction on the results obtained: max{/o, ^ {rj or Ja/S,}. If the defect is slightly split (for definiteness by 
magnetic field H) this condition turns into max{/o,|cL;|} ^ max{{?7 or Ja/^^}, gfisHS}. For nearly symmetrically 
coupled impurity one has: max{/o, |cli|} ^ max{{7] or Ja/^}, j^i — 52I}, where gi_2 are values of coupling with the 
corresponding sublattices (see Fig. ^| . 

The findings discussed above are vahd for isotropic interaction of the impurity spin S with AF: Hint = ffS(si -I- S2), 
where Si 2 are host spins from different sublattices. We consider also interaction containing only one component of 
S: Hint = gS^{sf + 82). The results in this case are quite specific. We show that xx- component of the impurity 
susceptibility is zero whereas yy- and zz- ones have only the non-resonant term. Our model is identical to the spin- 
boson one ^ without Hd if the interaction contains a term with only one component of S. The Hamiltonian can be 
diagonalized exactly-^ and an exact expression for xi'^) can be obtained. We perform the corresponding calculations 
and confirm the results obtained by our approach. One of the most interesting features of the exact result is that the 
static susceptibility has the form x(0) oc T~^~^, where C oc f^T/J. Within the first order of one has l/(4r)-term 
and the logarithmic correction. Thus, we see that in the modified spin-boson model taking into account the higher 
order logarithmic corrections results in the non-trivial power-law T-dependence of x(0)- 

The influence of the finite concentration n of the defects on the low-temperature properties of 2D AF is also 
considered. For not too small oj we find the logarithmic correction to the spin- wave velocity of the form nf^ In \J/u}\ 
and an anomalous damping of the spin- waves proportional to n/^jwl. Similar logarithmic correction to the velocity 
and damping were obtained in Rcf.^^«, where vacancies in 2D AF were studied. It is demonstrated that interaction 
of the spin waves with defects modifies the spectral function which acquires new terms proportional to n exhibiting 
weaker w-dependence. These terms should be taken into account at small enough lu and the problem should be solved 
self-consistently. The corresponding consideration is out of the scope of this paper. Within the range of validity of 
our study we do not obtain a renormalization of the magnetic specific heat which is proportional to in 2D AF 
without impurities. 

The results of the present paper can be applied to other systems with a degenerate defect in which the spectral 
fmiction is proportional to cj^ . 
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APPENDIX A: CALCULATION OF ImA^^(tj) IN 2D AF 

In this appendix we discuss properties of the imaginary part of the function A^i,(w) general expression for which 
is given by Eq. It is shown below that within the spin- wave approximation ImA^i,((jj) has the form © and 
expressions for the constant A, the characteristic energy Q and the tensor d^i, are obtained. 

We have from Eqs. ||3J) and iQ: 

^^^--H ^^J2^l + cos{kR,2Ms'l^s'Ch, (Al) 

k 

where {. . ■)uj denote retarded Green's function, TV is the number of spins in the lattice and R12 is the vector connecting 
two host spins coupled to the defect. Thus we have to calculate the spin Green's functions (s^,^sj^)tj of 2D AF. The 
Hamiltonian of 2D Heisenberg AF on the square lattice has the well-known form: 

H=jJ2^i^j- (A2) 

m 

We perform for beginning all calculations for isotropic Heisenberg 2D AF at T = and then consider the effect of 
finite T and of the additional small interaction stabilizing the long range order at finite T. 

Instead of dividing of the lattice onto two sublattices it is convenient to represent operators as follows (see, e.g., 
Refs»): 

Sk = xs^ + i/s^+ko + ^■^k+ko ' (A3) 



2 l«k 



(a^at^ 



2s 



is 



2s 



[a'a)k, 



(A4) 



where z axis is parallel to the magnetization of sublattices, x,y,z denote unit vectors directed along corresponding 
axes, ko = (tt, 0,7r) is the antiferromagnetic vector and s is the spin value. Substitution of Eqs. l|A3p and l|A4p to 
Eq. (|A2p leads to the following expression for the Hamiltonian: H = Eq + J2^=2^i^ where Eq is the ground state 
energy and Hi denote terms containing products of i operators a and a^. We consider in this paper the spin- wave 
approximation, i.e., we restrict ourself by the bilinear part of the Hamiltonian which has the form 



Ekal^au 



Bk 



'k'^-k 



flkfl-k 



(A5) 



where Ey^ = sJq, B]^ = sJ^ and = 2J(cosA;j. +cosfc2). As is seen from Eqs. (jAl|l . (|A3|I and (jA4|) . the only nonzero 
components of the spin Green's function are xx and yy and tensor has the form: 



Sf_„y, iffi,iy = x,y, 
0, if fi — z OT v : 



(A6) 



The corresponding components of A^t^(a;) can be derived using Green's functions 5(w,k) — {ak,al)uj, /(tj^,k) = 
(ak,a_k)w, 5(^7 k) = (aLk,a-k)w = 9*{-^,-^) and f''{uj,'k) = (a^^,ajj.)i^ = /*(-cj,-k). For two of them we have 
the Dyson equation: 



g{co, k) = g("\uj, k) + g(0)(L^, k)i?k/^(^, k), 
/t(u;,k)=g(0)(^,k)i?k5(c^,k), 



(A7) 
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where g''-^^ {uj,^) = (oj — £'k + iS) ^ is the bare Green's function. Solving Eq. ljA7|l one obtains: 

f^''^ =7 — r~IT2 2"' = — — — J, A8) 



where et = s -^Z Jq — is the spin- wave energy. As a resuh of direct calculations we have: 

r 1 

ImA^^(w) = -V— / dk [{1 + cos(kRi2)}(Jo - Jk) + {1 + cos((k + ko)Ri2)}( Jo + Jk)] — [S{uj - ek) - d{uj + eu)] , 

(A9) 

where the lattice constant is taken to be equal to unity and the integral is over the magnetic Brillouin zone. If 
^ sJ we have Jk « Jo(l — ek = cfc = VSsJk and cos(kRi2) « 1 — (kRi2)^/2. Notice that (koRi2) = tt 

mod 27r if the impurity is coupled to spins from different sublattices and both terms in the first square brackets in 
Eq. (|A9|) are proportional to . Then integration in Eq. (|A9|I can be easily carried out if one takes advantage of the 
approximation for magnons similar to Debye one for phonons: the spectrum is assumed to be linear, ek = cA;, up to 
cut-off momentum fee defined from the equation 2N = V(2tt)~^ J^"^ dkk, where V is the area of the lattice. As a 
result we lead to expression © for ImA^^(cj), where 

e = eke = SVt^sJ, A=—. (AlO) 

The factor A should be multiplied by 2 if the defect is coupled to four host spins (two by two from each sublattice). 

It is well known that there is no long range order in Heisenberg 2D AF at T > Oi^ Nevertheless it is shown 
theoreticallj*2iiLii and confirmed experimentalljii^ that the spin waves are well defined in paramagnetic phase of 2D 
AF if their wavelength is much smaller than the correlation length ^ cx exp(const/r). Thus, the above result for 
ImA^^(a') is valid when |w| 3> Ja/S,, where a is the lattice spacing. 

It is easy to conclude that if a small interplane interaction of the value of 77 ^ J is taken into account the above 
result for ImA^j,(cL') is valid when ^ 77 (see discussion in Sec. Ill A^ . At the same time ImA^^([x)) has another 
(jj-dependence if ^ rj. 

Finally, we note that when the impurity is coupled to one host spin we have A^y(aj) = Sk('^-k*k)'^ instead 
of Eq. (|A1|) . Comparing this equation with IjAip and ljA9|l one infers that the spectral function is proportional to a 
constant in this case. 



APPENDIX B: MATRIX STRUCTURE OF PSEUDOFERMION GREEN'S FUNCTION AND THE 

VERTEX FOR 2D AF 

In this appendix we discuss the matrix structure of the dressed pseudofermion Green's functions Gmm' (x) and the 
pseudofermion vertex Tpmm'ix + uj,x) for 2D AF. Some lower-order diagrams for the self-energy (w) and for 

^Pmm' are shown in Figs.Oland^l respectively. Let us discuss firstly the self-energy. Its matrix structure is determined 
by the corresponding products of operators S'^ and tensors d^^. For instance, this product for the second diagram in 
Fig. Ohas the form J^fifi'i^u' S^S'^S^ d^v'd^^i . One can make sure that such combinations are proportional to the 
unit matrix using the following evident representation of an arbitrary matrix A of the size 2x2 via Pauli matrices: 

A = ao + ((Ta), (Bl) 

where oq = Tr(A)/2 and a = Tr(crA)/2. It is seen from the view of the tensor d^^ given by Eq. ljA6p that the 
combinations of S"^ contain products of even number of matrices ax and <jy. According to HB1|I such combinations 
are proportional to the unit matrix. 

The similar consideration can be carried out for the vertex Tpmm'{x + uj,x), where P is one of Pauli matrices. 
The vertex matrix structure is determined by the products of operators 5"^, tensors d^^i, and P. For example, the 
corresponding product for the fourth diagram in Fig.^has the form S^S^PS^ S" d„^id^^i . It can be easily 

shown using ljBl|) that such combinations are proportional to P. 



APPENDIX C: CALCULATION OF THE IMPURITY SUSCEPTIBILITY 



We present in this appendix some details of the impurity dynamical susceptibility calculation. We use for this the 
general expression (|19|) and Eqs. (|26|l . H33|) and (|37|l for the Green's function and the branches of the vertex. The 
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expression for xp{^) is derived up to the order of In this order the interaction does not change the average number 
of pseudofermions Af given by Eq. IjlSI) . i.e. M = 2e~^^^ . To show this let us express Af as an integral of the Green's 
function»2£ 

1 f°° 2e^'^/^ 
Af^ V/ dxn{x)lmG^m{x) = / dxe-''/'^lmG{x) (CI) 

where n{x) = [e^l^ + 1)^^ is the Fermi function and G{x) is given by Eq. (|26|) . We make a shift by A in the last 
part of Eq. ICll) and replace n(x + A) by as it was done in Eq. ((T^ for xp('^)- Because 7(0;) and Z{x) are 

exponentially small at negative x \i\x\^ T (see Eq. (|2U) the integrand in Eq. I|C1|I does not increase exponentially 
as a; — > — cx). It is easy to make sure that the terms proportional to cancel each other in Eq. HC1|I . 

According to Eqs. (|19|l . and (|4()|l the dynamical susceptibility can be represented as a sum of three components. 
The first one, xi(ci;), originates from Eq. H19|) as a result of replacement of the vertex by unity. The second, xiiy^)^ 
appears from /-^-terms in Eqs. (|33() and (|40|l . The third, X3('^): is a result of replacement of the vertex by the third 
term from Eq. (|4()|l . 

The expression for Xil'^) can be brought to the form 

Xx(uj) = — / dxe-^/^[G(x + w) + G*(x-tj)]ImG(a;). (C2) 

The integrand in Eq. (|C2|I does not increase exponentially as x — s- — 00 because 7(0;) and Z{x) obey the property 
(I31|l . Using Eq. H26|l for Green's functions it is convenient to represent Eq. 1C2|I in the following form: 



Xi(co) = [Jo(cc') + Ji(co) + J2(co)] + [J*(-co) + Ji*(-C^) + J*(-L^)], (C3) 

2+^2' 

" -^7(2;) \Z{x) + Z* (x)] + x\Z(x) - Z* {x)\ ' 



M^) = r ^^^"^^^ 2^i"r2 . (C4) 



p2 /-oo 1 



47ri ,1^, X + uj + iff 







(C5) 



27r 7_oo a; + a; + .x^ + F^^ 

where in all denominators we replace j{x) and "f{x + co) by Iq, their values at \x\ T and |x + cij| ^ T, respectively 
(see Eq. (150)1 V It can be done because |7(a;)| <C \x\ at |a;| > T. 

The integration in Eq. (|C4I) can be easily carried out if one notes that the main contribution arises from the area 
of X ~ w, r, where w 1 — xjT (it will be clear soon that the second term in this expansion is essential). As a 

result we have: 

^o(")^-?^(l-^). (C7) 



2 w + 2iro V T 

To take the integral in Eq. (|C5(I for J\{lS) we consider a contour integral with the same integrand. The contour 
is presented in Fig. It consists of four lines which are parallel to the real axis. They pass through points x = 0, 
X = ilo — iJ, X = i/o + i(S and x = zttT. It can be shown using definitions of 7(x) and Z{x) that the integrand is an 
analytical function inside the contour. Note, the contour envelops the cut of Z(x)* passing through the point x = iJo 
along the real axis. As a result we have: 

r , , -^iT. ..J, , 7(a: + i7rr)(xi + X + iTTT) + 7(xi + x + i7rr)(x + mT) 
''^("^ = ---I^ / ^-'^-i^ ki|iV(xOA(xi) (.. + ..r)3(xi+x + z.r)^ 

(C8) 



27r ttG 7_oo (.X + j7rT)3(xi + X + zttT)^ 

P^/^%^-.ro/T r ^^^^^e-/^|xi|iV(xi)A(xi) r 11 11 



47ri 7r0 J-oo x + a; + 2i/o 



xi + X — «(5 X — i(5 xi + X + z(5 X + 2(5 



where the first term is the result of integration over the line of the contour passing through x = ittT and the second 
term describes the sum of integration over lines passing through x = iPo — iS and x = iFo + iS. It is seen from Eq. (|C8|) 
that the first term is of the order of and can be discarded. The second term in Eq. I|C8I) can be simply calculated 
using the equation (xi + x — iS)^^{x — iS)^^ — (xi + x + iS)^^{x + iS)^^ = 2TTi[S{xi + x)(x — iS)^^ + S{x){xi + iS)^^]. 
As a result we have: 

^ ' 2 uj + 2iro ttG x + Lo + 2iro ^ ' 
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One can carry out the integration in Eq. ljC6|l for J2(w) in the similar way. As a result we obtain that J2(u;) = Ji{lu) 
and Xi(cl;) has the form: 



Xi(^) 



p2 / 2zro 

2 \T{Lo + 2ira) 



1 



dx 



\x\N{x)k{x) 







sgn(a;)A(x) 
a; + w + 21/^0 



(CIO) 



(Cll) 



where we omit lo in the denominator of the integrand in the first term. 
The quantity X2('-^) can be expressed as follows: 

X2H = JH + j*(-c^), 

"p2 t2 j-oo 

J{tu) = T^-Ri'h^ / dxdxie~''/^\xi\xiN{xi)K{xi)G{xi+x + uj)G{x + uj)luY{G{x)G{xi+x)}. (C12) 
27r ttB 

The integral in Eq. lfnT2|) can be taken similar to those of Ji{ijj) and J2{i^), the result being 

X2(w) 



p2 

2 



2^ r j.mx)m p_ r- sgn(.)A(.) 



(C13) 



where we omit cj in the denominator of the integrand in the first term. 
The expression for Xai'-^) can be brought to the form: 



p2 



2 T{LU + 2ir) 



Zr 



dxe-'^'^Gix + uj)G*{x)K{x). 



(C14) 



The area near poles of the Green's functions is essential in this integral. Therefore, we can replace e ^/"^ by unity 
and K{x) by K{Q) = (/q — r)/^. As a result we have: 



p2 



2ir 



2iro 



P r , \x\N(x)K{x) p 
ttB a; + 2«/o 7r6 



_r(w + 2ir) T(w + 2iro) 

Summing Eqs. IjClOjl . ljC13(l and (|C15|I we lead to Eq. lEjl for the impurity susceptibility. 



dec 



|a;|A^(a;)A(x) 



X + 2z/o 



(C15) 



APPENDIX D: GREEN'S FUNCTIONS OF 2D AF WITH IMPURITIES 

We derive in this appendix Green's functions of operators a and considered in Appendix ^ in the case of 2D 
AF with impurities. To investigate the influence of the impurities we have to take into account the corresponding 
interaction Q in Dyson equations. It is assumed that N, Ni ^ oo so as Ni/N = n — const, where N and Ni are the 
number of spins in the lattice and the number of impurities, respectively. Within the linear spin-wave approximation 
the equations for g and /t have the form 

k) = .9(")(a>, k) - g^ns[l + cos(kRi2)][/(°)(^, k) + g^°\iu, k)]xA^M^: k) + p{iu, k)] 
+ g^ns[l - cos(kRi2)][/(°)(^, k) - g^°^{cj, k)]x^(c^)[5(^, k) - /t(L^, k)], 
Pico, k) - /t(0)(c^, k) - g'ns[l + cos(kRi2)] [.9<"' k) + /t(")(^, k)]xA^Moj, k) + p{uj, k)] 

+ g^ns[l - cos(kRi2)] [g^°^ (c., k) - /t(") (cj, k)]xy (^) k) - /t (c., k)] , (Dl) 

where superscript (0) denotes Green's functions without impurities which are given by Eq. I)A8|I and Xxi^^) and Xy{^) 
are the impurity susceptibility given by Eq. I|49|) with P — and S'^, respectively. Equations IDll) can be easily 
solved with the result 



5(w,k) = g(-cj,k)* = 
P{cj,k) = /(-c.,k)* = 



sJq + uj — 2g'^nsx±ito) 

-s Jk + 2g^nsx±{uj) cos(kRi2) 



where x±i^) = Xx{^) = Xy{^) and the denominator k) is given by Eq. 1)65(1 . 



(D2) 
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FIG. 1: Unit cells of 2D AF with (a) symmetrically and (b) asymmetrically coupled impurity spins are presented. Strengths 
of coupling with corresponding host spins g and g\ 7^ §2 are depicted. The local Neel order is also shown. Only symmetrically 
coupled impurities are discussed in this paper. 




^2 






FIG. 2: Lower-order diagrams for the impurity dynamical susceptibility xp(y^) a graphical representation of the result of 
the overall series summation. Lines with arrows represent the pseudofermion Green's functions. Wavy lines denote Green's 
functions of operators e^(Ro) of the host system (see Eq. (|5J). 
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FIG. 4: The zero-, the first- and the second-order diagrams for the pseudofermion vertex Vp(x -I- uj^x). 
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FIG. 3: The first- and the second-order diagrams for the self-energy E(a;) of the pseudofermion Green's function. 
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FIG. 5: Graphical representation of averaged interaction of one impurity with 2D AF. Bold solid and bold wavy lines de- 
note, respectively, dressed pseudofermion Green's functions and the boson Green's function renormalized by interaction with 
impurities. 




FIG. 6: Contour of integration used in Appendix |n| for the impurity susceptibility calculation. This contour consists of four 
lines which are parallel to the real axis. They pass through points x — 0, x = iFo — iS, x = iFo + iS and x — inT. 



